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ABSTRACT 

Waveform multigrid method is an efficient method for solving certain classes of time- 
dependent PDEs. This paper studies the relationship between this method and the analogous 
multigrid method for steady-state problems. Using a Fourier- Laplace analysis, practical 
convergence rate estimates of the waveform multigrid iterations are obtained. Experimental 
results show that the analysis yields accurate performance prediction. 
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1. Introduction. Waveform relaxation, also called dynamic iteration or Picard- 
Linderlof iteration [7], is a technique for solving ordinary differential systems of initial- 
value type [5] [8]. Its key idea is to solve an ordinary differential system in n variables 
by solving sequence of subsystems in fewer variables. This nature of the method allows 
independent integration with different timesteps for each of subsystems. Thus the 
method is usually considered in the context of parallel or multirate algorithms [2]. 

The convergence of the waveform relaxation methods has been analyzed by Miekkala 
and Nevanlinna [7]. They discussed this issue over infinite time interval and showed 
that, for linear heat equation with standard spatial discretization, the convergence rates 
of the waveform iterations are similar to those for the analogous steady-state problems. 
Therefore the convergence could be too slow for the waveform relaxation to be compet- 
itive with standard timestepping methods. 

Multigrid techniques (in space) can be incorporated to accelerate the convergence. 
This was studied by Lubich and Ostermann [6] who compared the multigrid performance 
of waveform iteration with that of static iteration (i.e., the iteration for steady-state 
problems) for one-dimensional heat equation. Combining the analysis on the smoothing 
rate of high frequencies, they conjectured the waveform multigrid performance for two- 
dimensional case. They showed that the typical multigrid acceleration can be achieved 
with an estimated convergence rate, which is similar but not quite as good as the one 
for the steady-state problems. 

A number of parallel waveform multigrid algorithms have been proposed and im- 
plemented [9] [11]. Analytical and experimental results have shown that the waveform 
multigrid methods can be implemented on a parallel computer with satisfactory effi- 
ciencies. 

In this paper we study the relationship between the waveform multigrid method for 
time-dependent PDEs and the standard multigrid for the corresponding steady-state 
problems. This study is important since the steady-state multigrid has beep investi- 
gated extensively while the properties of waveform multigrid algorithm are relatively 
unknown. We present a Fourier- Laplace analysis for obtaining practical convergence 
rate estimates of waveform multigrid iterations. The approachs used in this paper are 
simple, applicable to wide class of applications, and provide insight into the details of 
the basic interaction between the coarse grid correction and the waveform relaxation. 

This paper is organized as follows. Section 2 briefly introduces the waveform re- 
laxation method. In Section 3, the waveform relaxation and multigrid iteration are 
combined. A theorem is proved which indicates that the convergence rates of waveform 
multigrid are essentially the same as those for the analogous steady-state problems when 
number of smoothings is large. Section 4 gives the details of Fourier-Laplace analysis, 
which is used for obtaining exact convergence rates. As an example, analytical com- 
parison of multigrid convergence rates in practical number of smoothings is made for a 
two-dimensional heat equation. Finally, the comparison to the measured convergence 
rates is presented in Section 5. 

2. Waveform relaxation method. Waveform relaxation method was originally 
proposed for solving ordinary differential equations consisting of subproblems with few 
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external variables in circuit simulation [5]. Unlike standard time stepping methods, 
it iteratively partitions a large system into loosely coupled subsystems and integrates 
each subsystem independently. Hence, it is well suited for parallel machines, especially 
massive parallel machines. 

The method can be described as follows. Consider a linear initial value problem 


(1) — + Au = f, t > 0, u(0) = u 0 . 

Let A be split as A = M — N. Under certain conditions, Eq.(l) can be solved iteratively 
by 

(2) + Mu (v) = yv u (w_1) + /, t > 0, u (v) (0) = U 0 , 

tit 

which is equivalent to the integral equation 


(3) u ( ”> = Su< u - ,) + <£, 

where <5 is a linear integral operator on L P (R + , C n ) (1 < p < oo) with kernel k s (t) = e tM N 
and <f> is a function depending upon / and M. 

The general convergence results of the scheme (2) over entire interval t > 0 were 
given by [7]. In their work, the Laplace transformation was applied to the time variable 
t and it showed that the convergence rate for scheme (2) was given by the spectral 
radius of 5, derived as 


(4) p(S) = max p(S(z)), S(z ) : Laplace transform of S. 

We shall follow their approach and concentrate on the damped Jacobi relaxation and the 
red-black Gauss-Seidel relaxation for linear equations with time-independent coefficients 
of the form 


(5) 


du 

dt 


+ Lu = /, t > 0, u(0) = u 0 , 


where L is a linear elliptic operator. Let Lh be a discrete approximation of L and 
assume it can be written as 


( 6 ) 


L h =d 


I 

-B ' 

i 

i 

-R 

I 

• i 


Then the damped Jacobi relaxation and the red-black Gauss-Seidel relaxation have the 
Laplace transform 

d 

( 7 ) 


Sj{z) = {zI + Mj)-'Nj = 


z + d/u 


(i-1)/ B 
R (i-i)/ 


Mj = -/, Nj = Mj- L h , 

U> 


2 


0<u> < 1; 


Sgs(z) = [zl + Mgs) ' Ngs 


,_£_).[« (* 4 * ) B 

K z + d J 0 RB 


Mgs = d 


/ 0 

-R I 


Nc,s — Mas — Lfe- 


In the following sections, the subscripts of the matrices will be dropped when the context 
is clear. 

3. Waveform versus steady-state multigrid method. The convergence of the 
waveform relaxation method can be accelerated if the multigrid technique is incorpo- 
rated in space. In this section, we shall show that, the convergence rate of waveform 
multigrid iteration converges to that of steady-state multigrid iteration as the number 
of smoothings increases. In next section, a Fourier- Laplace analysis will be introduced 
to show that the performance of the waveform multigrid iteration for time-dependent 
PDEs is virtually as good &s the standard multigrid iteration for the corresponding 
steady-state problems for practical number of smoothings. 

The multigrid method is adapted to the waveform iteration for a time-dependent 
PDE in the following way. First, the equation is discretized in space to give a semi- 
discrete problem. Next, the multigrid iteration with waveform relaxation is applied 
to the space variables. As an example, a two grid V-cycle for Eq.(5) is illustrated as 
follows: 

• Perform v\ pre-smoothings: 

= Nu^ + /, t > 0, u< v >(0) = u 0 , 
at 

where Lh = M - N, v = 1,2, • • • ,Vj, u (0) = u (0) (<) is given. 

• Restrict the defect from grid h to grid H: 

( 9 ) d h := + L h u^ - /, d H :=lUh- 


• On the coarse grid, solve 


• Correct 


+ Lhv — du, v(0) = 0. 


= - I h H x 


where I is a suitable interpolation from grid H to grid h. 
• Perform V 2 post-smoothings on u. 
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Similar to the Full- Approximation-Scheme (FAS) formulation of multigrid, one can 
formulate a coarse grid problem as 

d u H r _ jHj , d ( , t rH,.M 

-^ + LhUh - -l h d h + —{l h u h ) + L,Hi h u h , 

u//(0) = /fuo 

and the correction step a 5 

« = “l”' 1 + - /fd 1 " 1 ) 

in order to handle non-linear problems. This however should be used in an FMG 
algorithm where the problems is solved first on coarse levels to obtain a good initial 
approximation to fine ones. 

The error = u *') — u of a complete two grid V-cycle iteration described above 
satisfies 

(10) e (<) = Ve { ‘ _1) 

for an integral operator V, which has the Laplace transform (see [6]) 

(11) V(z) = S(zr(I - l},(z + L H )- l l”(z + L h ))S(zy\ Rez > 0, 

where S(z ) is the Laplace transform of the smoother being used. In order to indicate 
the dependency of V and V(z) on the number of smoothings v = + v 2 , we use the 

notation V(u) and V(z,v) whenever it is necessary. 

Assuming that all the entries of V ( z ) are rational functions of z vanishing at infin- 
ity with poles having negative real part, and taking V as an operator on L P (R + ,C n ) 
(1 < p < cx>), the spectral radius p(V) = lim^oo ||V fc || 1/,fc satisfies ([6]) 

(12) p(V) = max o p(V(z)). 

Note that 5(0) and F(0) are respectively, the smoothing and multigrid two-grid cycle 
operators for the corresponding steady-state problem (or static problem) LhU — f. 

Theorem 1. The spectral radius for the waveform two-grid V-cycle operator 
V = V(u) for equation (5) satisfies 

(13) Jhn |p(V(u))-p(K(0,u))| = 0 

for either damped Jacobi or red-black Gauss-Seidel smoothings. 

Proof. We shall prove the case of the red-black Gauss-Seidel relaxation only. The 
proof for the damped Jacobi relaxation follows similarly. Let CG(z) denote the Laplace 
transform of coarse grid correction operator, 

CG(z) = I-I^iz + Ln^Ilfiz + L,) 

= CG(0) + CG a (z) 


( 14 ) 


with 


(15) CG^z) = zl k „(z + L H )-'(L- H 'I»U - /a"). 

Since the spectral radius 

p(S V2 (z)CG(z)S v '(z)) = p(CG(z)S v (z)), v = vj + v 2 , for all z, 
for simplicity, we rewrite the two-grid operator (11) as 


(16) 


V(z,v) = CG(z)S v (z) 

= CG{0)S v {z) + CG a (z)S v {z). 


As noted in Section 2, the red-black Gauss-Seidel smoother (8) 


(17) 


S'M = tiz) 


2v 


0 (^)5C ,V - 1 

0 C v 


C=RB, 


i 2v 


(_d_' 

l 2v 

' 0 

BC v ~ l 1 

V*+<L 

I 

0 

o 


Let z := zjd , we have 


V{z,v) 


(z + l) 


2v 


V(0,u) + A (z,v) 


with 

(18) 


A(z, v) = F^ v \z)A(z,v), 


F M (z) = 77 


(z+1) 2 *’ 


A(z,v) = CG(0) 


0 BC v ~ l 
0 0 




0 ( l + z)BC v ~ 1 

0 C v 


Given c > 0, there is a matrix norm || • || t , such that (see [4], pp.297) 


and 


Therefore 


p(V{z,v)) < \\ — ^2^ (0,^)^ + ||A( 2 ,U 


^(Q»^)lle < Pi jz + 


'(*+!) 


2v 


‘(*+0 


p(V(u)) < max — - • p(V(0, i>)) + e + max l|A(z, v) 

rV V ” ~ Rez>0 \z + l| 2tJ X V ” Rez">0 V ’ ’ 
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which leads to 


(19) 


0 < f>(V(.j)) - f>(V(0,v)) < ( + max IIAtz.^ll,. 


The matrix function A(z,v) is analytic in Rez > 0, continuous on Rez = 0, and 
A(oo,t>) is bounded for all v > 1 (because ^(O) — * 0 as v —* oo [7]), so does its e-norm 
||A(z,u)|| ( , thus 

(20) max ||A(z, v)|| ( < c(e), 


where c(e) is a function of e independent of z and v. Combining Eq.(18) and Eq.(20), 
we have 


max ||A(z, u)|| t < max 

Rez>0 K Rez> 0 




F {v) (±i 


1 


y/2v - 1 


c(e) 


—*• |F^(0)| • c(e) = 0, as v — * oo. □ 

4. Fourier-Laplace analysis. Fourier analysis has been used to provide exact 
convergence rate of multigrid iteration for some steady-state problems [1]. In order to 
analyze the convergence rate of waveform multigrid iteration, first, the Laplace trans- 
form is introduced in time variable to convert a m-dimensional semi-discrete linear 
differential equation with time-independent coefficients 

(21) u t + L h u = f 


u(0) = 0 

to equivalent algebraic equations defined over the right half complex plane 
(22) (z + L h )u{z ) = /(z), Rez > 0, 

Then multigrid method applied in space only is analyzed by considering its action on 
each of these equations. This is done using Fourier analysis for each equation in (22), 
and combining the results to obtain the convergence rate of waveform multigrid iteration 
for Eq.(21). In this section a Fourier-Laplace analysis is applied to the waveform two 
grid V-cycle iteration described in Section 3. 

Let us begin with the error formula (10) 

e« = Ve<K 


Its Laplace transform is 

(23) e (,) (z) = V(z)e (,-1) (z), Rez > 0. 
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Let 


9 = (0\O 2 ,---,O 2m ), 0 j = 0'(modir), j = 2, • • • 2 m , 
and X(0), V(0,z) denote the Fourier modes on grid h and the symbol of I/(z), i.e., 

X{9) = [exp(i9' x / h ) , • • • , exp(i$ 2m x / h)] 


V(0, z ) — [Vj,k(9,z)]2 m x 2 m 


satisfying 

(24) V(z)X(6) = X(0)V(9, z) for all 9. 

Then 

(25) p{V{z)) = sup/^V^z)) = max p(V(0,z)). 

0 \<2 

The symbol matrix V(9 , z) is of order 2 m , a tiny matrix comparing to V(z). Its spectral 
radius can be calculated accurately for each given 9 and z. Therefore />(V), the conver- 
gence rate of the waveform two grid iteration, is obtained by computing p(V(0, z) over 
0 < |0 T | < f and Rez — 0: 


p(V) = 


max max 
fle2=0 0<|e'|<a 


p( W*)b 


Example. Consider the heat equation on a square Cl = (0, 7r) x (0, tt) with Dirichlet 
boundary conditions 


u t — Au = /, ( t,x ) G (0, oo) x Cl 


u = g, (t, x) (E [0, oo) x dfl, u(0, x) = u 0 (x), x £ Cl. 
Let Lh correspond to the five point Laplacian 


(26) 




-1 


-1 

4 -1 
-1 


The Laplace transform of the two grid iteration operator V is given by (see Eq.(16)) 
(27) V(z ) = CG(z)S v (z) 


CG(z) = I - I h „(z + L h )-'!^z + L h ) 
7 


with 



and the corresponding smoother S'(z). Since the space spanned by the Fourier mode 
X{9) is invariant under each of operators in Eq.(27), we have 

(28) V{e,z) = (l-I h H {i)(z + L H (20))-'i«(0){z+Um^,z), Rez> 0, 


where ~’s indicate the matrix symbols and H = 2ft is assumed. For this particular 
example, the fine grid operator is represented as 


L k (0) 


U(9') 


U(9 4 ) \ 


with 


0 = (<Ma), 


the coarse grid operator as 


i (a \ _ 4(sm 2 (^i/2) + sin 2 {6 2 j 2)) 


L H (2e) = [L H ( 29 1 )}. 


The bi-linear interpolation is chosen for I fa and the restriction operator If? = {Ifa) T ■ 
Their matrix symbols are 






iH (9) = ( 


1 + cosO i x/ , 1 + cos0 2 x 
o X n /> 


i h H w = AW- 

Smoothers for the damped Jacobi and the red-black Gauss-Seidel relaxation are repre- 
sented as (see Eq.(7)-(8)) 


Sj(9,z) 


S(9\z) 


S(9 4 ,z) 




$<»,*) ~(d- u U(e)), 0 < w < 1; 

zu) + d h 


Sgs{ 9 , z) 


Sa( 9 ,z) 0 

o S b (0,z) 



1 + a. —(1 + a) 

1 — a —(1 — a) ’ 



1+6 -(1 + 6 ) 
1-6 — (1 - 6 ) J ’ 
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Table 1 

Damped Jacobi Waveform Mulligrid (v = 2/3) 


V 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

p{V(v)) 

.6626 

.4390 

.2909 

.1927 

.1486 

.1269 

.1124 

.1004 

.0908 

.0838 

rino,-)) 

.6626 

.4390 

.2909 

.1927 

.1332 

.1152 

.1003 

.0881 

.0793 

.0728 


Table 2 

Red-Black Gauss-Seidel Waveform Multigrid 


V 

i 

2 

3 

4 

5 

6 

7 

8 

9 

10 

p(V(v)) 

.2439 

.1567 

.1132 

.0861 

.0737 

.0614 

.0548 

.0473 

.0430 

.0390 

pcm*)) 

.2439 

.1517 

.1046 

.0765 

.0669 

.0543 

.0474 

.0404 

.0368 

.0336 


a = WT7j 7) {cose ' + CC39l) ’ 6 = + c ° s6l) ' 

Tables 1 and 2 present computed p(V(v)), the spectral radius of waveform multi- 
grid operator, and p(V(0,t;)) of corresponding steady-state multigrid operator for the 
damped Jacobi and the red-black Gauss-Seidel relaxation. As it shows, the estimated 
convergence rates of waveform multigrid iteration are virtually as good as those for 
steady-state problems. Since extensive research has been done in multigrid methods 
for steady-state problems [1] [3] [10], this analytic comparison is very useful in pre- 
dicting the performance of waveform multigrid methods. Although our discussion was 
done only for the model heat equation, one may expect the same performance of the 
waveform multigrid iteration as that of steady-state one in wide classes of applications. 
The treatment of general problems, e.g., non-constant coefficients or general domains, is 
done in a similar way using frozen coefficients argument which is applicable for smooth 
coefficient problems. Its theoretical rigorous justification is involved and needs the use 
of pseudo-difference calculus. 

5. Comparison and Conclusion. We have shown that the estimated conver- 
gence rate of waveform multigrid iteration is almost undistinquishable from that of 
steady-state one. Now, we compare the experimental results of the convergence rate to 
the analytic ones discussed in last section. 

Consider the two-dimensional heat equation on Q = (0, ?r) x (0,7r) 

tt< — A u = /, (t,x) G (0, tf] x 


u(i,x) = 0, (t,x) G [0,t/] x dft, 


u(0, x) = u 0 (x), x G fi, 

where u 0 was generated randomly to excite all possible Fourier modes. The space 
derivatives were discretized by central differences with uniform fine grid size h x . For 
time integration, as in [7], trapezoidal rule was used. The step size was always chosen 
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Table 3 

Damped Jacobi Waveform Multigrid 


V 

analytic p(V(v)) 

measured ^(V(i>)) 

worst case 

avareage case 

1 

.6626 

.6657 

.6403 

2 

.4390 

.4401 

.4160 

3 

.2909 

.2964 

.2693 

4 

.1927 

.1890 

.1836 


Table 4 

Red- Black Gauss-Seidel Waveform Multigrid 


V 

analytic p(V(v)) 

measured /?(V(t;)) 

WM^-l) + 0(hl) 

worst case 

avareage case 

1 

.2439 

.2560 

.2282 

.2500+0(/^) 

2 

.1567 

.1137 

.0740 

.1625+0(/i2) 

3 

.1132 

.0761 

.0524 

.1295+0(/»2) 

4 

.0861 

.0670 

.0437 

.11 lO-l-O(A^) 


as h t = .01 on [0, tf]. Note, the efficiency of time integration is not the concern of this 
paper. 

The experiments were run for two-grid V-cycle with v damped Jacobi (cj = 2/3) 
or red-black Gauss-Seidel relaxations. The fine grid mesh size in space was h x = tt/ti, 
n = 8,16,32,64. Although the spectral radius p(V(i>)) was studied on entire time 
interval t G [0,oo), finite interval [0,f/] had to be used in experiments. For each set 
of tests, we used both tf = 1 and tf = 10. However, we found that the measured 
convergence rates did not depend on the size of time interval. 

To measure the spectral radius p(V(u)), we used the asymptotic ratio of the defects 
(see Eq.(9)) 

maxi Il4 ,+1) ll2 

max, \\4% 

Using the matrix split of L^, the calculation of the derivatives in the defects can be 
avoided. Because of extensive computations involved, ||dj ^||2 was evaluated only at 
t = tf since the convergence of the waveform iteration is determined by the error at 
t = tf. 

Tables 3 and 4 present the comparison results. They show that the analytic con- 
vergence rates obtained via Fourier- Laplace analysis almost coincide with the real ones, 
a result evidenced for the steady-state problems [1], As a reference, Table 4 also lists 
values of \\frjol 2v — 1), the bounds of p(V(v)) proved for one-dimensional problem by 
Lubich and Ostermann [6]. They conjectured the bounds for that of two-dimensional 
heat equation as |^A/o(2w — 0 + Observing our analysis in Section 4, with 

standard discretization of L, nine-point restriction and bi-linear interpolation, /9(V(u)) 
is independent of the grid size h x . The major difference of our approach to the one 
in [6] is that, instead of using eigenvectors of L we used Fourier modes, which are 
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much easier to manipulate, extendible to high-dimensional problems and wide class of 
applications. Most important of all, this approach is able to give the exact convergence 
rates for special model problems and sharp estimates for general problems. 
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